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Abstract 

Standard quantum state reconstruction techniques indicate that a detection efficiency of 0.5 is an absolute 
threshold below which quantum interferences cannot be measured. However, alternative statistical techniques 
suggest that this threshold can be overcome at the price of increasing the statistics used for the reconstruction. 
In the following we present numerical experiments proving that quantum interferences can be measured 
even with a detection efficiency smaller than 0.5. At the same time we provide a guideline for handling the 
tomographic reconstruction of quantum states based on homodyne data collected by low efficiency detectors. 


1 Introduction 


Homodyne detection is an experimental method that is used to reconstruct quantum states of coherent light by 
repeatedly measuring a discrete set of field quadratures [TJ|^. Usually, a very high detection efficiency and ad-hoc 
designed apparatuses with low electronic noise are required [^. New methods capable of discriminating between 
different quantum states of light, even with low detection efficiencies, will pave the road to the application of 
quantum homodyne detection for studying different physical systems embedded in a high noise environment 
For this purpose, specific quantum statistical methods, based on minimax and adaptive estimation of the Wigner 
function, have been developed in |&11 . These approaches allow for the efficient reconstruction of the Wigner 
function under any noise condition, at the price of acquiring larger amounts of data. Hence, they overcome 
the limits of more conventional pattern function quantum tomography |12f[19| . The important consequence 
of this novel statistical approach is that the 0.5 detection efficiency threshold can be overcome and quantum 
tomography is still practicable when the signals are measured with appropriate statistics. The scope of this paper 
is to report the results of this method tested by performing numerical experiments. Indeed, we consider a linear 
superposition of two coherent states and numerically generate homodyne data according to the corresponding 
probability distribution distorted by an independent Gaussian noise simulating efficiencies lower than 0.5. By 
properly expanding the set of numerically generated data, we are able to reconstruct the Wigner function of 
the linear superposition within errors that are compatible with the theoretical bounds. Our results support the 
theoretical indications that homodyne reconstruction of linear superposition of quantum states is indeed possible 
also at efficiencies lower than 0.5. 


2 Wigner function reconstruction 

Let us consider a quantum system with one degree of freedom described by the Hilbert space (K) of square 
integrable functions ijj{x) over the real line. The most general states of such a system are density matrices p, 
namely convex combinations of projectors onto normalised vector states 

p = Xj \'4’j){'4’j\ ) Aj > 0 , Aj = 1 . 
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Any density matrix p can be completely characterised by the associated Wigner function Wp{q,p) on the 
phase-space {q,p) G namely, by the non positive-definite (pseudo) distribution defined by 


Wp{q,p) = 


{2nf 


dMdz;e*(“«+’'P)Tr 


pe 
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= — due*“P {q-v/2\p\q + v/2) 
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( 1 ) 


Here q and p are the position and momentum operators obeying the commutation relations [q, p] = i, h = 1, 
and \q ± v/2) are eigenstates of q: q\q ± v/2) = {q± v/2)\q ± v/2). Notice that Wp(q,p} is a square integrable 
function: 

( 2 ) 


27r f dqdp \Wp{q,p)\'^ = Tr (p^) < 1 . 
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Among the advantages of such a representation, is the possibility of expressing the mean value of any operator O 
with respect to a state p as a pseudo-expectation with respect to Wp{q,p) of an associated function 0{q,p) over 
the phase-space, where 

( 3 ) 
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Indeed, by direct inspection one finds 
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In homodyne detection, a monochromatic signal photon state is mixed with a coherent reference state, a so-called 
local oscillator, by a 50/50 beam splitter. The output is collected by two photodiodes and the difference 
photocurrent is measured. It can be proved that, when the local oscillator is s ign ificantly more intense than 
the signal, the homodyne photocurrent is proportional to the signal quadrature 20 . Denoting by d and d^ the 
single mode annihilation and creation operators associated with the signal, the quadrature operator is defined as 

de~^^ + 


^(f) — 


V2 


( 5 ) 


where is the relative phase between signal and local oscillator. The continuum set of quadratures with (/ G [0, tt] 
provides a complete characterization of the signal state. Using the annihilation and creation operators d, d^ one 
constructs position and momentum-like operators, q = {d + a^)/-\/2 and p = {d — d^)/{i'/2). With respect to the 
latter, the quadrature operator reads: 


= q cos (j) + p sin (j) . (6) 

Quadrature operators have continuous spectrum extending over the whole real line, Xfj)\x) = x \x); given a generic 
one-mode photon state associated with a density matrix p, its diagonal elements with respect to the (pseudo) 
eigenvectors 

Pp{x,(t)) := {x\p\x) , (7) 

represent the probability distribution over the quadrature spectrum. 

In homodyne detection experiments the collected data consist of n pairs of quadrature amplitudes and 
phases (A^, ^i)-. these can be considered as independent, identically distributed stochastic variables. Given the 
probability density Pp{x, </), one could reconstruct the Wigner function by substituting the integration with a 
sum over the pairs for a sufficiently large number of data. However, the measured values x are typically not the 
eigenvalues of x^, rather those of 


+ 



y , 0 < ?7 < 1 , 


( 8 ) 


where y is a normally distributed random variable describing the possible noise that may affect the homodyne 
detection data and p parametrizes the detection efficiency that increases from 0 to 100% with p increasing from 
0 to 1 1^. The noise can safely be considered Gaussian and independent from the statistical properties of the 
quantum state, that is, y can be considered as independent from x^p. Then, as briefly summarised in Appendix 
A, the Wigner function is reconstructed from a given set of n measured homodyne pairs {Xf, ^e), ^ = 1, 2,..., n, 
by means of an estimator of the form 


= ^lniy^P)Xr{q,P) , ^IJq^p) = ^ ^ (^[( 9 , P); 4*^] “ 


( 9 ) 
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This expression is an approximation of the Wigner function in (27) by a sum over n homodyne pairs (Xi, 

The parameter h serves to control the divergent factor exp ( 7 ^^), while r, through the characteristic function 
Xr(,q,p) of a circle C'r(O) of radius r around the origin, restricts the reconstruction to the points {q,p) such 
that . Both parameters have to be chosen in order to minimise the reconstruction error which is 

conveniently measured |10| by the L^-distance between the true Wigner function and the reconstructed one, 
||Wp — 1 T^’^|| 2 - Since such a distance is a function of the data through W^'^, the L^-norm has to be averaged 
over different sets, M, of quadrature data: 


K:niP)=E \\W,-Wl 


h,n 


= E 


dqdp 


/B2 


Wpiq,p) - 


( 11 ) 


where E denotes the average over the M data samples, each sample consisting of n quadrature pairs {X(, 
corresponding to measured values of with </> S [ 0 , 7 r]. In |10| , an optimal dependence of the parameters r and 
h upon the number of data, n, is obtained by minimizing an upper bound to r 


3 Interfering Coherent States 

Homodyne reconstruction is particularly useful to expose quantum interference effects that typically spoil 
positivity of the Wigner function: it is exactly these effects that are claimed not to be accessible by homodyne 
reconstruction in presence of efficiency lower than 50%, namely when 77 in (|^ is smaller than 1/2 |15| . However, 
in it is theoretically shown that 77 < 1/2 only requires increasingly larger data sets for achieving small 
reconstruction errors. However, this claim was not put to test in those studies as the values of 77 in the considered 
numerical experiments were close to 1 . 

Instead, we here consider values 77 < 1/2 and reconstruct the Wigner function of the following superposition 
of coherent states 

I'I'a) = : , |a) = e-l“l'/2e““' |0) , (12) 

y'2(l + e-2|“P) 

with a any complex number ai + m 2 ^ C. 

The Wigner function corresponding to the pure state pa — |^a)(^a| is shown in Figurefor a € M. Its 
general expression together with that of its Fourier transform and of the probability distributions Pp{x^ (/>) and 
p'^{x^(j)) are given in Appendix B. 
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Figure 1: Wigner function corresponding to the pure state pa = |^a)(^a| (o^i = S/v^; 0^2 = 0). 

In Appendix C, a derivation is provided of the L^-errors and of the optimal dependence of h and r on the 
number of data n and on a parameter P that takes into account the fast decay of both the Wigner function and 
its Fourier transform for large values of their arguments. The following upper bound to the mean square error in 
0 is derived: 

< A , A = 4 Ai(7) + , (13) 

nh 

^ T he functional relation between the parameters h and r on n also depends on an auxiliary parameter 0 > 0. This was introduced 
in [l0| to characterise the localisation properties on of the Fourier transforms of the Wigner functions of the following class of 

density matrices addressed in that context: A^^s,l = 'Ip • J 1^ [tFp] (g,p)p ^ < (27r)^L|-. 
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with 0 < /3 < 1/4 and 7 as given in (26). 

As explained in Appendix C, the quantities Ai 2,3 do not depend on h, r and n. By taking the derivatives 
with respect to r and h, one finds that the upper bound to the mean square deviation is minimised, for large n, 
by choosing 


1 

''~h~ 


I logn 


(14) 


We generated M = 10 samples of n = 16 x 10® quadrature data {X(, $£) distributed according to the noisy 
probability density p^p{x,(j)) explicitly given in (32) of Appendix B, considering an efficiency lower than 50% 
(77 = 0.45). Starting from each set of simulated quadrature data we reconstructed the associated Wigner function 


by means of ([^ and (10). The averaged reconstructed Wigner functions E 
in Figure]^ for two different values of the parameter /3. 


h,7 


{q,p) for r] = 0.45 are shown 




^= 0.05 



Figure 2: Averaged reconstructed Wigner functions E over M = 10 samples of n = 16 x 10® noisy 

quadrature data (efficiency ry = 0.45). Two different values of /3 are considered. 


The mean square error of the reconstructed Wigner functions has been computed as in (11) and compared 
with the mathematically predicted upper bounds A. The dependence of the upper bound reconstruction error 
on the parameter f3 is discussed at the end in Appendix C. In Table we compare the reconstruction errors 
A/’(((p) with their upper bound A for two significant values of /?. 


p 

<] 

A 

0.05 

0.081 

2.39 

0.1 

0.076 

26.07 


Table 1: Calculated A/’0(p) for M = 10 samples of noisy quadrature data (77 = 0.45) for two different values of 
p. Comparison with the mathematical prediction of the upper bound A. 


Despite common belief, the interference features clearly appear in the reconstructed Wigner function also for 
efficiencies lower than 50% and the reconstruction errors are compatible with the theoretical predictions. In the 
next section, we make a quantitative study of the visibility of the interference effects. 
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3.1 A witness of interference terms 

The interference effects in the state jl'a) can be witnessed by an observable Oa of the form 

A _ |a)(-a| + I - a)(a| 

2(l+e-2|«l=) ■ 

With respect to an incoherent mixture of the two coherent states, 

Pa\ = A |a)(a| + (1 - A) I - a){-a\ , 0 < A < 1 , 

its mean value is given by 

"^T:\Pa\Oaj = ^_2la\2 ' 

Therefore, from (§ it follows that the phase-space function Oa{q,p) associated to is 


Oa{q,p) = 


e cos ( 2 \/ 2 (qa 2 +pai)) 

■ v^(l+e-2|“l=) " 


Q. — Q-i -|- io^2 • 


(15) 

(16) 

(17) 

(18) 


For details see (28) and in Appendix B. 

Let us denote by lF“^g^(( 7 ,p), the estimated Wigner function Wjl’^{q,p) in (|^ for the j-th set of collected 
quadrature data. It yields a reconstructed mean value 


{Oa)j,rec= / dq dp O^iq,p) p) 


(19) 


of which one can compute mean, Av(< Oa >rec), and standard deviation, Sd(< Oa >rec), with respect to the 
M sets of collected data: 


Av{{6a)rec) = 


Sd{{Oa)rec) = 
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( 21 ) 
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We computed Av{{Oa)rec) and Sd(< Oa >rec) with M = 10 simulated sets of noisy data with p = 0.45 for two 
different numbers of simulated quadrature data (see Figure]^. We repeated the procedure for different values of 
the parameter /3. The results are presented in Figure]^ where the error bars represent the computed Sd((OQ,)rec)- 

In order to be compatible with the interference term present in iTa), the reconstructed Wigner functions 
should yield an average incompatible with the incoherent mean value in (|17[), namely such that 


Av(< Oa ^rec) 


,-2|a| = 


l + e-2|“P 


A Sd(<C Oa ^rec) ■ 


( 22 ) 


We thus see that the condition in (221 is verified for p = 0.45, that is the reconstructed Wigner functions are not 


compatible with incoherent superpositions of coherent states, if enough data are considered. We also notice that 
the same behavior is valid for the high efficiency p = 0.95. 

The dependence of the errors on /? can be understood as follows: when /3 decreases the integration interval 


in (101 becomes larger and approaches the exact interval [—c»,-|-oo]. Nevertheless this occurs at the price of 
increasing the reconstruction error. This can be noted both in Figure(larger error bars) and in Figure]^ 
(increasingly noisy effects in the reconstructed Wigner function). This problem can be overcome with a larger 
number of data samples M, that reduce the reconstruction noise and compensate for the effect of decreasing /3. 


4 Conclusions 

We simulated quadrature data corresponding to high electronic noise and detection efficiencies lower than 0.5. 
Under these operating conditions the Wigner function of a linear superposition of two coherent states could 
be reconstructed using the tomographic techniques developed in j^. Moreover, by taking into account the 
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Figure 3: Av{{Oa)rec) — ^ function of j3. The error bars represent Sd ( ( Oa)rec )• For each /3, 

M = 10 set of n noisy quadrature data have been considered. The square markers refer to rj = 0.45 (n = 16 x 10® 
blue marker and n = 5 x 10® green markers) while the round ones refer to rj = 0.95 (n = 16 x 10®). The error 
bars for 77 = 0.95 have been multiplied by 20 in order to make them more visible. 


decay properties of the Wigner function along with those of its Fourier transform, we have checked that the 
numerical reconstruction errors are compatible with the theoretical error bounds computed there. Furthermore, 
the reconstruction of the quantum interference pattern of the Wigner function has been supported by a numerical 
study of the variance of an operatorial interference witness excluding that the reconstructed interferences might 
be an artefact of the reconstruction algorithm. 

We have thus confirmed 1) that, as theoretically predicted in |9 10 , a 0.5 detection efficiency is not, as 
often stated in the quantum optics literature, an absolute threshold below which homodyne quantum state 
reconstruction is generically impossible and 2 ) that, instead, by suitably enlarging the size of the set of collected 
quadrature data, and using alternative techniques different from standard pattern function quantum tomography, 
one may indeed access quantum features even in low efhciency conditions. 

These results also provide the tools for quantum state reconstruction in these operating conditions and set 
the boundaries of the applicability of the novel statistical approach to homodyne quantum state reconstruction 
by checking the increase of the number of quadrature data necessary for faithful reconstruction with decreasing 
detector efhciency. 
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A Wigner function reconstruction 


The quadrature probability distribution Q can be conveniently related to the Wigner function by passing to 
polar coordinates u = ^coscp, v = ^ sincj), such that 0 < (j) < tt and —00 < ^ < + 00 : 
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where F\pp{x, ^)](0 denotes the Fourier transform with respect to x of the probability distribution: 


/ +00 

da:e"“«Pp(a;,(?i) . 
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( 24 ) 


Since y can be considered a normally distributed random variable independent of the noise affected distribution 
of the eigenvalues of in ([^ is given by the following convolution: 
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Its Fourier transform is connected with that of Pp(x,4>) according to 

F[pp(x,((i)](C) =e^«V[p^(x,(/))](^/v^) , with 7 := 


1 - V 
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By inserting F[pp{x into (241, one can finally write the Wigner function in terms of the noisy probability 

distribution p^{x, (p): 

p-K p + OO \p\ p-\-00 

Wpiq,p)= dp dxe*««“^^+P-"^-"/v^)e^«%p(x,(^) . (27) 

Jo J— OO J— OO 

B Coherent state superposition: Wigner function 

The Wigner function corresponding to the pure state pa = |4'q,)(4('q| and its Fourier transform read 
w (n rd — __ 7p-(iJ-%/2ai)^-(p-V2a2)^ i p-(g+\/2Qi)^-(p+%/2a2)^ 
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(29) 


For a generic Wigner function Wp{q,p) one computes the quadrature probability density Pp(x, p) in ([^ by means 
of the so-called Radon transform: 


(x|/3|x) = / dp Wp(xcos(/) —psin(/), Xsine/)-f pcos0) , 
Jr 


(30) 


It follows that the probability density Pp[x,4>) and the noise-affected probability density p)J(x,(/)) in (251 are 
given by: 


Pa(x,(^) = 
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a{4)) = tti cos (j) + a 2 sin p , /?(</') = 0:2 cos p — ai sin p . 
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C Upper bound reconstruction error estimation 

Here we derive an upper bound to the mean square error of the reconstructed Wigner function. This analysis 
is necessary in order to find an optimal functional relation between the free parameters in the reconstruction 


algorithm such to minimise the reconstruction error. For this purpose we follow the techniques developed in 10 
adapting them to the case on a linear superposition of coherent states. Using (|^ and (101, one starts by rewriting 
the error in (111 as the sum of three contributions: 
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C°(0) denoting the region outside the circle C'r.(O), where . The first and the third term correspond 

to the variance and bias of the reconstructed Wigner function respectively, while the second term is the error 
due to restricting the reconstruction to the circle (7^(0). 

Given a density matrix p, the second term can be directly calculated. This is true also of the bias; indeed, 
because of the hypothesis that the pairs (Xe,^e) are independent identically distributed stochastic variables, it 
turns out that 


E 


WIM^P) 


1 




kI[UpW-^ 


= ^E 

TT 


A'd [(«.!>); *1-^) 


pTZ \f\ /•+00 

/ Act) / d^T;^ / 

0 -J—Xlh J — oo 


(36) 


differs from the true Wigner function Wp{q,p) in (271 by the integration over ^ being restricted to the interval 
[—1/h, 1/h]. Moreover, its Fourier transform reads 
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where X[-i/h,i/K\{\\w\\) is the characteristic function of the interval [—1/h, l/h]. Then, by means of Plancherel 
equality, one gets 
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The variance contribution can be estimated as follows: firstly, by using ([^ and (101, one recasts it as 
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Then, a direct computation of the first contribution yields the upper bound 
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with o(l) denoting a quantity which vanishes as h when /i —^ 0. On the other hand, the second contribution can 
be estimated by extending the integration over the whole plane {q,p) S and using (371 together with (|^: 
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Let us consider now the specific case oi p = pa, the superposition of coherent states defined in (12). The 
auxiliary parameter (3 labelling the class of density matrices Ap^s.L in footnote p] with s = 2 can be used to 


further optimize the reconstruction error A^’^(/3 q,). In particular, since X] 7 =i j\ < M we get the 

upper bounds 
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which simultaneously hold for 0 < /3 < 1/4. 

Then, by means of the Cauchy-Schwartz inequality, one can estimate the contribution (|34| to the error. 
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Altogether, the previous estimates provide the following upper bound to the mean square error in (33l-(35): 

AlliM < A , A = 4 Ai(7) + A,iP) , 


(51) 


where Ai^ 2,3 do not depend on h, r and n and Ai( 7 ) = •y7r/(167r^y^) is the leading order term in (40). By 
setting the derivatives with respect to r and h of the right hand side equal to 0, one finds 
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Whenever /3 is such that the arguments of the logarithms are much smaller than the number of data n, to leading 
order in n the upper bound to the mean square deviation is minimised by 


1 / log n 

h y /3 + 27 ■ 


(54) 


The range of possible values of /3 is 0 < /3 < 1/4. However, the upper bound A becomes loose when /? —>■ 1/4 and 
/3 —> 0. In the first case, it is the quantity A 3 (/ 3 ) which diverges, in the second one, it is the variance contribution 
which diverges as the logarithm of the number of data. It thus follows that the range of values /3 G [/3q,/3i] where 
the numerical errors A/’(^(/3 q-) are comparable with their upper bounds A is roughly between /3o = 0.04 and 
Pi = 0.10 for rj = 0.45 as indicated by the following Figure]^ 





Figure 4: Upper bound reconstruction error A as a function of the parameter /3. Two efficiencies r] are considered. 
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